Accurate Calculation of Green Functions on the d-dimensional Hypercubic Lattice 
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We write the Green function of the d-dimensional hypercubic lattice in a piecewise form covering 
the entire real frequency axis. Each piece is a single integral involving modified Bessel functions of 
the first and second kinds. The smoothness of the integrand allows both real and imaginary parts 
of the Green function to be computed quickly and accurately for any dimension d and any real 
frequency, and the computational time scales only linearly with d. 
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Lattice Green functions arise in numerous problems in 
combinatorics, statistical mechanics, and condensed mat- 
ter physics. In the language of condensed matter physics, 
the Green function G(r, uj) is the propagator for a parti- 
cle of energy (or frequency) uj to move a distance r in a 
tight-binding model on a lattice with a nearest-neighbor 
hopping. Of particular interest is the local Green func- 
tion G(uj) = G(r = 0,ui), whose imaginary part is the 
local density of states. Closed-form expressions exist for 
G(uj) on the chain, square, cubic, honeycomb, triangular, 
body-centred cubic (bcc), face-centred cubic, diamond, 
and d-dimensional hyper-bcc lattices in terms of elliptic 
integrals and hypergeometric functions (see Ref. [T] fo r a 
useful review), as well as on certain fractal lattices)^ but 
in general it is necessary to perform integrals or infinite 
sums. 

In this work we address the problem of computing the 
local Green function of a d-dimensional hypercubic lattice 
with nearest-neighbour hopping t = , j 



G d (Lu) = 



1 



dk\ . . . dkd 



ui — (cos k\ 



cos k d ) 



i0+ 
(1) 



(where we have included an infinitesimal shift in the de- 
nominator as is conventional in condensed matter the- 
ory). This function has a branch cut along the real 
axis for ui G [— d, d], representing a continuous spectrum 
(band) of particle excitations with bandwidth 2d. 

There are several ways of computing explicit values 
for Gd(ui). The density of states of the d-dimensional 
hypercubic lattice, Ad(uj) = — ^lmGd(uj), is the convo- 
lution of Ad ± and Ad 2 where d = d% + d?. Since A d is 
known in closed form for d = 1, 2, 3, this allows Ad to be 
computed as a single convolution for d ~ 4, 5, 6, but for 
d > 7 two or more nested integrations are required. It 
is simple to derive the La urent series for Gd(uj) from a 
multinomial expansion)^ but due to the branch points 
at ±d, the Laurent series only converges for |w| > d, on 
or outside a circle of radius d (barring techniques such as 
Borel summation). Gd(co) may be written as a Fourier 
integral involving Bessel functions (in which case the ac- 
curacy of numeri cal integration is limited by oscillations 
of the integrand) or as an Laplace transform integral 
involving the modified Bessel function Iq (which only 
converges for frequencies outside the band edge) P This 



work presents integral representations for both real and 
imaginary parts of G d {ui) applicable on the entire real 
axis, both inside and outside the band; the integrands 
are smooth and well-behaved, which makes it easy to ob- 
tain numerical results to high precision. 
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FIG. 1. (a) Structure of the Green function in the complex 
plane, illustrated for d — 3. Thick lines represent branch 
cuts and dots represent branch points. The Laurent series 
converges only on and outside the dashed circle. In contrast, 
the method of this paper applies for any lu on the real axis, 
(b) Structure of Eq. Q in the complex plane. The branch cut 
runs below the negative real axis, terminating in a logarithmic 
branch point. The integration path runs along the positive 
real axis, and is deformed to run along the positive or negative 
imaginary axis depending on the value of u> + d — 2n. 

Equation ([T]) gives the Green function as a mul- 
tidimensional integral. Fourier-transforming to the 
time-domain, performing the wavevector integrals, and 
Fourier-transforming back to the frequency domain leads 
naturally to a single-integral formula in terms of Bessel 
J functions: 



G d (t) 



dui 
71 dh 



l G d (w) 



■ dk d 



6(i)e 



— it (cos fci + ...+cos kd) 



= -i9(t)J (t) d , 



G d {uj) = 



dt 



(2) 



From Eq. ^ one can quickly estimate G £ j(w) on the real 
axis using ordinary quadrature or fast Fourier transform 



2 



methods!*^ However, the oscillatory nature of the inte- 
grand limits the accuracy of this approach to a few digits. 
In situations like this, a useful strategy is to choose an 
integration path within the complex plane along which 
the integrand is smooth and non-oscillating. 6 However, 
for frequencies within the band (\lu\ < d) this trick is not 
directly applicable because the integrand grows exponen- 
tially in both the upper half-plane (UHP) and lower half- 
plane (LHP). Therefore, we split the Bessel function into 
Hankel functions H(t) = H { a 1] (t) and h(t) = H^ 2) (t), 



G d (uj) = -i J dt e , (3) 

and perform a binomial expansion to obtain 

■ d „oo 

^H = -pE(n)/ dte^H(t) d - n h(tr. (4) 

n=0 •'° 

The asymptotic behavior of the integrand, ignoring 
power-law factors, is e^^+rf- 2 ™). if w + d - 2n > 0, the 
integrand decays exponentially into the UHP, so Jordan's 
lemma allows us to deform the integration contour into 
the positive imaginary axis. Conversely, if uj + d— 2n < 0, 
we can deform the contour into the negative imaginary 
axis. Thus, 



G d (u) 



d 

- Vf c 



2'i 

0{uj + d-2n) 
+ @(-{u) + d-2n)) 



lut H{t) d - n h{i) n 
dt e lu}t H(t) d - n h(t) n 



Substituting t — it and t = —ir respectively into the 
two integrals collapses both integrations onto r € (0, oo), 
allowing us to combine integrands: 



The integrand in Eq. ^ has integrable logarithmic sin- 
gularities at r — 0. When uj e {— d, —d + 2, . . . , d}, corre- 
sponding to van Hove singularities, the integrand exhibits 
power-law decay at large r; otherwise, it decays exponen- 
tially. These asymptotic behaviors are easily dealt with 
using standard transformations such as those built into 
Mathematical integration routines. At intermediate r 
the integrand is a smooth function that can be integrated 
quickly and accurately. 

For computational purposes, we can further optimize 
Eq. |5]) by explicitly writing out the real and imaginary 
parts of the Hankel functions. For r > 0, 



H(ir) = -iK(r), 
H(-zt) = -iK(t)+I(t), 

h{ir) = iK(r)+I(T), 
h(-ir) = iK(r), 



where for convenience we have defined K(t) = ^Kq(t) 
and J(r) = 2I (t) where K and Iq are modified Bessel 
functions (as implemented in Mathematica) . Substitut- 
ing in these relations, performing a binomial expansion, 
and rearranging summations leads to an integral contain- 
ing a sum of d+ 1 products of modified Bessel functions, 



G d (w) 



1 



dr 



^ E C 3m K d ' m {T)I m {T) 

m=0 
d-j-1 

e " T E D ]m K d - m {T)I m {r) 



(6) 



1 f°° _ _ 

1 J n=0 

Q(u + d - 2n) e - ur H{iT) d - n h{iT) n 
- 0(-(w + d - 2n))e" T H{-iT) d - n h(-iT) r 



The sum can be split into two terms, one for n < j and 
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one for n > j, where — L ^^^2^ J i s a staircase function: 



where the coefficients 



Gjm E („) ( m 

n—m 

d-j-1 

- E (n)( 
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n\jd-\-m—2n 



(7) 
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GdH = ¥ J Q dr 



J2 OH{irf- n h{irT 



n=0 



^ E OH{^r) d - n h(^rT 

n=j+l 



(5) 



can be precomputed. It should be noted that the above 
sums can be evaluated in closed form using the properties 
of binomial coefficients. As an illustration, setting d = 3 
leads to an integral formula for the cubic lattice Green 
function, 
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FIG. 2. Real and imaginary parts of hypercubic lattice Green functions for d — 1,2, 3, 4, 5, 6, 7 (shifted by — O.lrf for clarity). 
The van Hove singularities are almost unnoticeable for d > 5. As d — > oo, ImG(w) tends to a Gaussian as expected. 



G 3 (u) 



r-J 3 e 



dr < 



3~LJT 



u < -3 



-31 + iK(2K 2 coshwr - M 2 e UT ) -3 < u < -1 
sinhwr — 4iif 3 cosh wt 



3IK 2 e~ 

J3 e -cur 



iK(2K 2 coshwr - 3Pe~ UT ) 



-1 < w < 1 
1 < u < 3 
w > 3 



(8) 



(where r arguments of K and J have been suppressed 
for clarity). Eq. ^ agrees with existing formulas for 
M > 3 pi and it can be implemented with good ac- 
curacy. For example, G 3 (0) = -0.8964407887768..., 
which agrees with the closed-form resullP to 12 digits. 
The Green functions for d = 1 up to d — 7 are plotted in 
Fig. @ 

An elegant feature of this method is that it gives a 
piecewise representation of the Green function, where van 
Hove singularities emerge naturally from sudden changes 
in the integration path. This is in contrast with Lau- 
rent, Fourier, and Chebyshev series methods, where the 
summation proceeds blindly and requires a large num- 
ber of terms to approximate the singular behavior (if it 



converges at all). 

The approach can easily be generalized to off-diagonal 
Green functions G(r, oS). However, because it relies on 
the factorization property embodied in Eq. ([2]), there is 
no obvious way to generalize it to non-hypercubic lat- 
tices. 

Calculating properties of lattices in d > 3 dimensions 
is admittedly an academic pursuit with rather contrived 
applications (for example, Gq(uj) is the local two-particle 
spectral function of a cubic lattice). Nevertheless, the 
work presented here is valuable as a pedagogical example 
as how complex analysis can vastly improve the accuracy 
of numerical integration. The insights obtained from this 
exercise may well be useful in more practical areas. 
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